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Abstract. We study the frustrated Heisenberg model on the bilayer honeycomb lattice. The 
ground-state energy and spin gap are calculated, using different bosonic representations at 
mean field level and numerical calculations, to explore different sectors of the phase diagram. In 
particular we make use of a bond operator formalism and series expansion calculations to study 
the extent of dimer inter-layer phase. On the other hand we use the Schwinger boson method 
and exact diagonalization on small systems to analyze the evolution of on-layer phases. In this 
case we specifically observe a phase that presents a spin gap and short range Neel correlations 
that survives even in the presence of non-zero next-nearest-neighbor interaction and inter-layer 
coupling. 


1. Introduction 

The study of the possible disordered ground states on two-dimensional antiferronragnets has 
received a great interest in the last years. In particular, the existence of quantum disordered 
phases has been studied in the phase diagram of antiferromagnets in a single layer honeycomb 
lattice m m m m 0 m m m m ma. However, the study of the influence of possible interlayer 
coupling on these phases is scarce mmmm- From the experimental side, a significative 
progress on the study of the bismuth oxynitrate, Bi 3 Mn 40 i 2 (N 03 ), has been made by Smirnova 
et al. [15]. In this material the Mn 4+ ions form a honeycomb lattice and two layers of such 
honeycomb lattices are separated by bismuth atoms, forming a bilayer structure. The study 
of the magnetic susceptibility indicates two-dimensional magnetism and no long-range ordering 
down to 0.4 K, suggesting a nonmagnetic ground state[H>l [TB]. In addition, density functional 
studies indicate that dominant interactions are the interlayer interaction Jj_ and the nearest- 
neighbor interaction Ji on each layer [17]. 



Figure 1. (Color online) Schematic representation of the coupling interactions in 
Bi 3 Mn 40 i 2 (N 03 ). Colored areas correspond to the unit cells. The sites in each unit cell are 
labeled as (1,^4), (2, A), (1 , B) and (2 ,B). 


The aim of this paper is to study the zero temperature ground state of the frustrated spin- 
1/2 Heisenberg model on the bilayer honeycomb lattice. We study the 5 = 1/2 case where 
the quantum fluctuations becomes more important in order to characterize the quantum phases 
in the model. On the other hand, although the material Bi 3 Mn 40 i 2 (N 03 ) has 5 = 3/2, the 
substitution of Mn 4+ in Bi 3 Mn 40 i 2 (N 0 s) by V 4+ may lead to the realization of the 5 = 1/2 
Heisenberg model on the honeycomb lattice. We use two different mean held self-consistent 
approaches based on bosonic representations of the spin operators to study this system, combined 
with Lanczos and series expansion methods to support the mean held results. 

2. Self consistent calculations on the bilayer Model 

We study the following Heisenberg model on the bilayer honeycomb lattice 

H = J aAry)Sa(r) ■ Spi?) (1) 

r,r',a,(3 

where, S a (r) is the spin operator on site a corresponding to the unit cell r. a takes the values 
a = (1, A), (2, A), (1, B), (2, B ) corresponding to the four sites on each unit cell as depicted in 
Fig- HI together with the couplings J a fi{f,r*). The coupling J± does not introduce frustration 
in the system and then, at the classical level and T = 0, it does not affect the classical Neel 
phase, present for J 2 /J 1 < 1/6. In the quantum case the situation is much subtle, increasing 
J |_ Neel order is likely to melt giving rise to a non-magnetic phase. 

For large values of J± we expect the ground state to be an interlayer valence bond crystal 
(IVBC) with corresponding spins from both layers forming dimers (as ilustrated in Figure [T]). 
This limit represents an excellent starting point for the bond operators formalism and series 
expansion calculations. 

On the other hand, starting from the magnetically ordered phase, the Neel order can be 
destroyed both by increasing the frustration on each layer or increasing the coupling between 
layers. The destruction of Neel order in a frustrated single layer honeycomb lattice has been 
studied by means of various approaches mmmmmmmmm- In the following we use 
two different bosonic representations of the spin operators to study the influence of the interlayer 
coupling in the ground state of the bilayer model. 

2.1. Bond operators Mean field approach. 

First, we use the well known bond-operator method to study the bilayer antiferromagnet 
described by Hamiltonian ©• We start by introducing a bond-operator representation of spin 
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where i] = 1,2 is the layer index, S“^(r) is the spin operator on the sublattice A of layer r/ 

corresponding to the unit cell located at r (See figured]). Operators s^(r) and aj,(f) create 
singlet and triplets states (out of a vacuum |0)) in the vertical bond placed in sub-lattice A and 
are defined as: s}jO> = ^=(| Tl) - | |t», a£|0> = -^(| tt) - I U», 4l°) = ^(1 it) + I U», 
a||0) = -^(| tl) + | lt))> and similar expressions for sublattice B. These kind of representations 
where proposed by Sachdev]24] in order to treat quantum phase transitions between Neel and 
dimerized phases. Operators belonging to the same unit cell satisfy the bosonic commutation 
relations whereas operators belonging to different unit cells commute. The restriction that 
the physical states are either singlets or triplets leads to the constraints in each unit cell, 
aUa = 1 and + Y^ a b » b « = 1- Introducing the bond-operators representation 
of the spin operators in dT]) we obtain a bosonic version of the Hamiltonian. We transform 
Fourier and retain terms up to second order to write H = Hj_ + Hi + H^, where 

H_l = ~J±s 2 N+ ^^{a} x (k)a a (k) + bl x (k)b a (k)} , ( 4 ) 

k,OL 

Hi = yp y~] {tW (a a (k)b a (-k) + a a (fc)bt,(fc) + b^(fc)a a (fc) + (5) 

k,a 

+ 7 (k) (h a (k)a a (-k) + a^(fc)b a (k) + a^(fc)b^(-fc) + b a (£)at(£)) j , 

Ha = (2s 2 -5)iVA + A^{aL(fe)at(fe) + bL(fe)4(fe)}, (6) 

k,a 


where, we have assumed that condensation of singlets occurs, i.e. (s Q (r)) = s, -y(k) = 
s 2 (l + e lk ' ei + e lk ' e2 ) : e*i and e ?2 are the primitive vectors on a triangular lattice and A is a 
Lagrange multiplier related to the constraint in the number of bosons. Diagonalization by using 
a Bogoliubov transformation allows us to write the following expression for the ground state 
energy 


f = (2» 2 - 5)A - A/±(2s 2 + 1) + Y, f TT (‘A’’® + t '“ S>4) ) ’ <?) 

where u^\k) and u^\k) are triplet energies. The parameters s 2 and A are determined by 
solving self-consistently the saddle point conditions = 0 and = 0 which are used to 
evaluate the energy and gap of the system and compare with numerical techniques. 

In order to complement our study, we have performed series expansion (SE) calculations, 
starting from the limit of isolated dimers connecting spins from both layers via Jj_. Notice that, 
this kind of expansions remain valid in the same limit that the bond operator approach (i.e. in 
the limit of strong inter-layer coupling.) To this end we have performed a continuous unitary 
transformation on the original Hamiltonian, using the flow equation method. This technique 
allows to obtain perturbatively an effective Hamiltonian that keeps the block diagonal structure 
of decoupled dimers. We refer for details of the method to ref. p5]. 
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Figure 2. (Color online) Ground state energy per dimer as a function of J\ obtained by 
the self consistent bond-operator approach (blue circles), Lanczos (ED) on a 24 sites system 
(red squares) and 0(5) Series Expansion (SE) (yellow rhombi). Inset: triplet gap (same set of 
parameters as main panel) BO (blue circles) ED (red squares) and SE (yellow rhombi). 


For the present model we have performed 0(5) and 0(4) SE in J\ 2 for ground state energy 
and for triplet dispersion, respectively. Explicit expressions are too long to be printed explicitly 
but are available electronically upon request. In Figj2]we show the ground state energy per site 
as a function of J\ and J 2 = 0, obtained by BO (blue circles), 0(5) SE (yellow rhombi) and 
ED (red squares) on a system of 24 sites. As it can be observed, all the techniques predict an 
energy decreasing with the coupling of interlayer-dimer via J\. Furthermore, there is an excellent 
quantitative agreement between the three methods for small values of J\. On the other hand, 
triplet gap is shown in the inset of Figi2] for the same set of parameters as the ground state 
energy. Here we also observe that all the techniques predict a tendency to a closure of the gap, 
when Ji is turned on. Our calculations shows that BO, ED and SE predict the same behavior. 

2.2. Self consistent Schwinger Boson Mean-Field Theory 

As we have seen previously, bond operator and series expansion methods are both suitable 
to study the interlayer-dimer phase. In order to investigate the evolution of on-layer phases 
as a function of inter-layer coupling we apply a representation of the spin operators in terms 
of Schwinger bosons [57], S a (r) = ^ba(r) • a • b a (r). Here b a (r)^ = (b^(r),b^(f)) is a 
bosonic spinor corresponding to the site a in the unit cell at position r, a are Pauli matrices, 
and the constraint in the number of bosons Yla btuo-^b^^ = 2 S has to be satisfied on 
each site. In order to perform a mean field decomposition, we define the following SU{ 2) 
invariants, A a p{x,y) = \ ob a ^(x)hp _ a (y) and B a p(x,y) = \ bi i<7 (x)bp t - l7 (y). This 
decomposition allow us to treat ferromagnetism and antiferromagnetism on equal footing and 
has been successfully used to describe a number of quantum frustrated antiferromagnets[1] [2, 
is na m m mi. We perform a Hartree-Fock decoupling where the mean field parameters 
are the expectation values of the SU(2) invariants A and B. The mean field equations for 
these parameters A a p and B a p and the constraints in the number of bosons must be solved 
self-consistently (see refs. in, 0 and M for further details). To obtain the phase boundary 
between the magnetically ordered and disordered phases using the self consistent Schwinger 











Figure 3. (Color online) a) Phase diagram in the J 2 — J± plane obtained with SBMFT. 
Gray region corresponds to the Neel phase whereas green region corresponds to magnetically 
disordered phases, b) Phase diagram of the single layer case corresponding to Ref. [3j. Inset: 
Z3 order parameter corresponding to the line J 2 = 0.38 


boson mean field theory we study the boson spectrum. In the gapless region the excitation 
spectrum is zero at k = 0 , where the boson condensation occurs, this is characteristic of the 
Neel ordered phase. On the other hand, in the gapped region, the absence of Bose condensation 
indicates that the ground state is magnetically disordered. 

In Fig. [3j-a) we show the phase diagram in the J 2 — J± plane. For Jj_ J 2 one can 

expect a IVBC ground state adiabatically connected with the limit of decoupled dimers, i.e. 
two singlets per unit cell, between spins and 2,A(2,B) (see Fig. [[]). In the region 

0.2075 < J 2 < 0.289 there is a reentrant effect. In this range, Neel phase separates from J 2 axis, 
leaving a tiny space for a magnetically disordered phase. In this way, Neel phase is here not 
only limited by some value Jj_(J 2 ) from above, but also by a second value Jj)*(J 2 ) from below 
(See Figure 0. 

On the other hand, in the range 0.3732 < J 2 < 0.398 ( J± = 0), there is evidence of the 
existence of an on-layer valence bond phase [3] (see Figure [3]-b). In this phase, SU(2) and 
translational symmetries are preserved, but Z 3 symmetry is broken. By turning on J± the system 
moves to the IVBC where the Z 3 symmetry is recovered. In the inset of Figure [3]-a, we depict 
the Z 3 directional symmetry-breaking order parameter p (defined in [T 8 ]) vs J±. The behavior of 
this parameter suggest a Z 3 symmetry restoring. Finally, in the region 0.289 < J 2 /J 1 < 0.3732 
the ground state preserves SU(2), lattice translational and Z 3 symmetries and the spin-spin 
correlations are short ranged j!4] . This agrees with the evidence of a spin liquid phase in the 
phase diagram corresponding to J± = 0 M. 




















3. Conclusions 

In summary, we have studied the phase diagram corresponding to a frustrated Heisenberg model 
on the bilayer honeycomb lattice, by means of bosonic mean field theories, complemented 
with exact diagonalization and series expansion, and described the behavior of the quantum 
phases as the interlayer coupling is increased. Using the Schwinger boson description we have 
determined the region where the system is Neel ordered and the lines where the Neel order 
is destroyed. We have determined an intermediate region where the phase diagram shows 
signatures of a reentrant behavior and we observe that for values of the interlayer coupling 
between (0.289 < J 2 /J 1 < 0.3732) the Neel order is absent at J± = 0 and the system presents a 
nonzero spin gap, whereas in the region (0.3732 < J 2 /J 1 < 0.398) each layer presents a nematic 
disordered phase[3]. In all the range of J 2 studied, the system presents signatures of an interlayer- 
valence bond crystal (IVBC) phase for J±/J\ > 4. This phase evolves adiabatically from the 
limit of decoupled interlayer-dimers. This is corroborated by bond operators self-consistent 
calculations and series expansions starting explicitly from the limit of isolated interlayer dimers. 
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